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1.  SUMMARY 


Knowledge  of  the  intensity  of  electromagnetic  radiation  in  the  radiation  belts  is  a  key  component 
to  understanding  the  lifetimes  and  overall  behavior  of  the  energetic  particles  existing  there. 

These  particles  can  damage  spacecraft  as  well  as  human  tissue  exposed  to  them.  Very  Low 
Frequency  (VLF)  radiation  from  lightning  is  thought  to  be  an  important  contributor  to  this 
interaction.  Here  we  attempt  to  construct  a  climatology  of  VLF  wave  energy  from  lightning  in 
low  earth  orbit,  which  can  be  used  to  estimate  its  intensity  throughout  the  radiation  belt  region. 
The  final  distribution's  spatial  and  temporal  variations  are  compared  with  the  in  situ  data  and 
found  to  capture  all  of  the  salient  features.  This  represents  the  first  global  climatology  of  this 
type  and  should  be  used  in  situations  where  estimates  of  the  VLF  intensity  from  lightning  are  a 
necessary  component  of  a  mission  or  analysis. 

2.  INTRODUCTION 

It  is  important  to  understand  the  intensities  and  origins  of  VLF  radiation  in  the  inner 
magnetosphere  primarily  in  order  to  characterize  the  sources  and  sinks  of  energetic  electrons  in 
the  radiation  belts.  The  radiation  belts  contribute  to  reducing  satellite  lifetime  in  orbit,  as  well  as 
being  a  primary  hazard  to  human  health  during  space  travel.  Precipitating  electrons  from  the 
radiation  belts  are  thought  to  affect  the  chemistry  of  the  middle  atmosphere,  e.g.  [1,2],  and  may 
contribute  to  the  poorly  understood  sensitivity  of  the  earth’s  climate  to  solar  variability  as 
exemplified  by  the  Milankovitch  cycles.  Previous  work  constructing  climatologies  of  VLF 
radiation  from  satellite  data  have  shown  their  utility  in  enhancing  our  understanding  of  these 
important  sources  and  sinks,  e.g.  [3,  4,  5].  Constructing  a  climatology  of  VLF  from  lightning  can 
help  us  to  categorize  and  separate  terrestrial/non-terrestrial  and  anthropogenic/non-anthropogenic 
sources  and  therefore  contributions  to  radiation  belt  particle  lifetimes.  In  addition,  understanding 
the  contribution  of  lightning  to  VLF  intensity  in  low  earth  orbit  (LEO)  is  important  in  order  to 
help  separate  and  characterize  the  global  propagation  characteristics  of  VLF  radiation  both 
within  the  earth-ionosphere  waveguide  and  through  the  ionosphere. 

The  magnetosphere  represents  the  area  near  the  earth  whose  composition  is  not  immediately 
dominated  by  the  flow  of  the  solar  wind,  but  is  partially  shielded  by  the  deflecting  effects  of  the 
earth’s  magnetic  field.  Within  the  inner  magnetosphere  the  Earth's  magnetic  field  dominates  and 
can  be  reasonably  represented  without  an  external  field  model.  This  region  generally  extends  to 
around  6  Earth  radii  (Re~6372  km).  Lying  within  the  inner  magnetosphere  are  the  radiation 
belts,  where  the  motion  of  energetic  electrons  (~1  keV  -  5  MeV)  are  confined  by  the  earth’s 
magnetic  field.  The  radiation  belts  are  roughly  toroidal  in  shape,  following  lines  of  force  of  the 
geomagnetic  field.  The  inner  and  outer  radiation  belts  are  separated  by  a  'slot'  region  where  there 
are  fewer  energetic  particles.  The  location  of  this  slot  is  energy  dependent  and  has  provided  a  test 
for  theories  attempting  to  explain  the  existence  of  the  radiation  belts.  The  inner  radiation  belt  is 
relatively  stable  and  transient  injections  have  been  observed  to  occur  only  during  strong 
geomagnetic  storms,  e.g.  [6,  7].  Energetic  electron  densities  in  the  outer  radiation  belt  are 
relatively  dynamic  and  more  closely  linked  to  solar  activity.  The  electron  slot  occurs  around  L  = 
2.5  -  3.5  and  depends  on  electron  energy  in  such  a  way  as  to  suggest  that  loss  rates  due  to  wave- 
particle  interactions  are  dominant  within  that  region  [8].  The  inner  radial  edges  of  the  radiation 
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belts  are  regions  where  the  atmospheric  neutral  density  is  large  enough  to  make  Coulomb 
interactions  the  dominant  sink  for  energetic  electrons. 

VLF  radiation  plays  a  central  role  in  the  plasma  dynamics  of  the  inner  magnetosphere.  The 
energetic  electron  distribution  function  in  the  region  is  generally  thought  to  be  controlled  by 
wave  particle  interactions  and  Coulomb  scattering.  Doppler  shifted  cyclotron  resonant 
interactions  between  whistler  mode  plasma- waves  (10s  of  Hz  up  to  10s  of  kHz)  and  energetic 
electrons  in  the  radiation  belts  can  cause  pitch  angle  scattering  in  the  particles.  This  can  lead  to 
deeper  penetration  towards  the  earth  along  their  bounce  paths,  which  leads  to  enhanced  Coulomb 
scattering.  Coulomb  scattering  causes  energy  loss  which  in  turn  enhances  the  rate  of  Coulomb 
scattering.  This  process  drives  a  cascade  ending  in  the  deposition  of  all  of  the  energetic  electron's 
energy  in  the  upper  atmosphere;  a  process  known  as  electron  precipitation.  Wave-particle 
interactions  can  also  lead  to  the  growth  of  whistler  waves  in  a  process  known  as  amplification 
due  to  resonant  instability  [9].  Despite  decades  of  study  it  has  proven  difficult  to  show 
conclusively  what  the  dominant  source  and  loss  processes  are  for  energetic  electrons  in  various 
parts  of  the  magnetosphere,  in  part  because  of  a  lack  of  understanding  about  the  sources  and 
sinks  of  the  waves  involved. 

Theoretical  treatments  of  energetic  electron  lifetimes  within  the  radiation  belts  rely  on  empirical 
estimates  of  the  electron  distribution  functions  and  lifetimes,  but  also  on  estimates  of  the  wave 
energy  thought  to  exist  there.  In  seminal  work  by  Able  and  Thome  [10,  11,  12],  estimates  were 
made  for  the  wave  intensity,  occurrence  rate,  and  propagation  angle  of  all  waves  then  thought  to 
contribute  to  particle  lifetimes.  These  included  Plasmaspheric  Hiss  (estimated  from  work  by 
Thorne  et  al.,  [13];  and  Tsurutani,  Smith,  and  Thorne,  [14]),  Lightning-Generated  Whistlers 
(estimated  from  Burgess  and  Inan,  [15];  and  Draganov  et  al.,  [16]),  and  VLF  Transmitter  Signals 
(estimated  from  Inan,  Chang,  and  Helliwell,  [17]).  Utilizing  bounce-averaged  quasi-linear  pitch 
angle  diffusion  coefficients,  the  equilibrium  distribution  functions  and  precipitation  lifetimes 
calculated  were  found  to  accurately  describe  characteristic  structures  observed  in  the  radiation 
belts.  This  conclusion  was  recently  called  into  question  by  a  reanalysis  of  the  transmitter  signal 
estimates  by  Starks  et  al.  [18,19],  which  suggested  that  median  field  strengths  were 
overestimated  by  roughly  20  dB  and  so  could  not  impact  electron  lifetimes  as  necessitated  by  the 
Able  and  Thorne  analysis.  Due  to  the  coupled  nature  of  these  calculations  this  suggests  that  the 
other  assumptions  made  by  Able  and  Thorne  [10]  should  be  revisited. 

Lightning  generated  whistlers  (LGW)  begin  as  radio  frequency  electromagnetic  impulses  created 
by  the  currents  in  lightning  strokes.  These  impulses  are  commonly  known  as  sferics.  Sferics  are 
converted  into  plasma  waves  at  the  lower  ionospheric  boundary  and  from  there  obey  the 
conventions  of  a  plasma  wave  in  a  magnetized  cold  plasma.  That  is,  the  propagating  modes  are 
right  hand  circularly  polarized  waves  below  the  local  electron  cyclotron  and  plasma  frequencies, 
and  above  the  ion  cyclotron  frequency.  They  travel  along  trajectories  detennined  by  the  local 
magnetic  field  and  plasma  properties  (e.g.  electron  density  gradients  and  positive  ion 
populations).  They  are  known  as  whistlers  because  their  frequency  range  overlaps  with  the  audio 
and  they  can  sometimes  be  detected  by  radio  receivers  on  the  ground  after  having  traveled 
through  the  inner  magnetosphere.  On  this  journey  they  become  dispersed  in  frequency  and  so 
'whistle'  characteristically,  with  higher  frequencies  arriving  first  [20]. 
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In  order  to  be  observed  from  a  ground  station  it  is  thought  that  whistlers  must  be  guided  by  a 
field  aligned  density  irregularity  or  ’duct'.  This  is  because  Ducts  constrain  the  motion  and  wave 
nonnal  of  the  VLF  radiation  to  be  along  the  field  line,  and  the  propagation  angle  with  respect  to 
the  external  magnetic  field  near  a  strong  gradient  in  refractive  index  (like  the  top  or  bottom  of 
the  ionosphere)  must  be  near  zero  or  reflection  will  occur.  Non-ducted  waves  become  oblique 
over  time,  that  is  their  wave  normal  becomes  transverse  to  the  local  magnetic  field,  so  they  do 
not  penetrate  the  ionosphere.  Non-ducted  VLF  radiation  undergoes  spectral  dispersion  in  space 
as  well  as  time,  as  the  guiding  effect  of  the  magnetic  field  is  frequency  dependent.  Many  non- 
ducted  whistler  trajectories  include  frequency  dependent  inversions  to  their  group  velocity 
direction  with  respect  to  the  geomagnetic  field  known  as  magneto  spheric  reflections.  These  are 
distinct  from  reflections  occurring  at  strong  gradients  in  the  refractive  index  and  rely  on  the 
closing  of  the  refractive  index  surface  for  a  wave  with  frequency  near  the  local  lower  hybrid 
surface  [21].  Whistlers  that  undergo  these  events  are  known  as  magnetospherically  reflected 
(MR)  whistlers.  A  given  whistler  may  undergo  a  number  of  reflections  before  dropping  below 
detectability. 

Whistlers  are  often  characterized  based  on  how  many  reflections  they  have  experienced.  A 
whistler  heading  directly  up  through  the  ionosphere  to  a  satellite  location  is  called  a  0+,  or 
fractional  hop,  whistler  and  one  headed  down  at  a  satellite  (or  ground  station)  after  passing  the 
magnetic  equator  once  being  called  a  1  -  (or  single  hop)  whistler.  Both  ducted  and  nonducted 
whistlers  have  been  observed  in  the  magnetosphere  [22]  but  space  based  observations  of  ducted 
whistlers  are  rare.  The  fact  that  multi-hop  whistlers  can  be  observed  from  the  ground  testifies  to 
the  stability  of  some  field  aligned  ducts.  The  issue  of  what  fraction  of  whistler  energy  is  ducted  is 
central  to  calculating  their  importance  to  radiation  belt  electron  lifetimes.  This  is  because  the 
phase  direction  of  ducted  whistlers  remains  field  aligned,  while  that  of  non-ducted  whistlers 
becomes  oblique  with  respect  to  the  local  geomagnetic  field  and  though  their  group  velocity  is 
guided  by  the  field  direction  they  do  not  follow  it  closely.  Oblique  whistlers  are  individually  less 
efficient  at  producing  the  pitch  angle  scattering  which  results  in  electron  precipitation  [23],  but 
these  may  be  more  numerous  or  have  a  small  effect  in  a  wider  range  of  situations  [24].  The 
relative  prevalence  and  climatology  of  ducts  in  the  plasmasphere  and  inner  magnetosphere  is  still 
an  active  area  of  research. 

Lightning-induced  electron  precipitation  (LEP)  events  associated  with  whistlers  have  long  been 
detected  from  rockets  [25],  satellites  [26],  and  from  the  ground  [27].  Observing  correlations 
between  satellite  based  measurements  of  precipitating  electrons  and  ducted  whistler  observations 
at  Palmer  Station,  Antarctica,  Voss  et  al.  [28]  suggested  that  detectable  LEP  events  are  driven 
primarily  by  ducted  whistlers,  supporting  previous  work  by  Burgess  and  Inan  [15]  which 
analyzed  only  ground  based  observations.  Lauben,  Inan,  and  Bell  [24]  and  later  Bortnik,  Inan, 
and  Bell  [29]  focused  on  nonducted  whistler  wave  energy  and  concluded  that  it  may  play  a  role 
in  the  formation  and  maintenance  of  the  slot  region.  Inan  et  al.  [30]  interpreted  satellite  data  on 
precipitating  electron  flux  timing,  number,  and  energy  in  the  context  of  a  nonducted  interaction 
model  and  showed  reasonable  agreement.  That  work  also  observed  large  regions  associated  with 
thunderstorms  containing  enhanced  precipitating  electron  flux,  and  suggested  that  they  are 
produced  and  maintained  by  lightning.  These  regions  showed  up  as  large  scale  enhancements  to 
the  overall  background  electron  density,  and  not  discrete  impulsive  events.  These  observations 
emphasize  that  it  is  important  to  understand  the  overall  climatology  of  LGW  and  how  they 


3 


couple  into,  and  are  transported  within,  the  magnetosphere;  as  well  as  from  the  perspective  of 
individual  LEP  events.  This  broader  perspective  was  also  shown  to  be  useful  in  work  by 
Gemelos  et  al.  [31],  which  demonstrated  a  seasonality  in  electron  precipitation  consistent  with 
lightning  and  inferred  that  lightning  is  a  significant  contributor  to  the  formation  of  the  slot 
region. 

Using  satellite  observations  to  infer  that  specific  precipitation  events  are  driven  by  specific 
observed  whistlers  has  proven  difficult  because  sferics  in  the  earth  ionosphere  waveguide  can 
generate  whistlers  out  to  at  least  1000  km  [32,33],  and  even  that  appears  to  be  an  observational 
rather  than  theoretical  limit.  Sferics  are  routinely  observed  on  the  ground  at  distances  of  6000  km 
[34]  or  more  from  their  originating  lightning  strike  depending  on  the  frequency  range 
investigated  [35].  It  similarly  becomes  more  difficult  to  positively  identify  a  sferic  observed  on 
the  ground  with  a  whistler  observed  in  orbit  as  the  distance  between  the  observations  grows. 
Propagation  of  the  sferics  within  the  earth  ionosphere  waveguide  is  nonunifonn,  as  is 
transionospheric  absorption.  Additional  low  altitude  and/or  columnar  ionization  drives  additional 
adsorption;  situations  favored  during  the  daytime  and  at  low  latitudes.  It  also  seems  that  there  are 
locations  where  the  transformation  from  a  sferic  to  a  whistler  is  more  or  less  efficient  such  that 
observation  at  certain  locations  is  more  likely,  e.g.  [36].  A  final  complication  occurs  because  the 
electromagnetic  pulse  generated  by  a  lightning  strike  can  actually  modify  the  ionosphere  as  it 
passes,  generating  optical  emissions  such  as  'elves',  causing  local  heating  of  electrons,  and  even 
driving  ionization  which  in  turn  changes  the  local  refractive  index. 

3.  METHODS,  ASSUMPTIONS,  AND  PROCEDURES 
2.1  Organizational  Framework 

The  intent  of  this  study  is  to  construct  a  climatology  of  whistler  wave  energy  from  lightning  at 
660  km.  The  approach  followed  here  is  to  start  from  previously  constructed  climatologies  of 
lightning  near  the  earth’s  surface.  That  distribution  is  then  considered  as  a  source  of  VLF  energy, 
which  can  propagate  within  the  earth-ionosphere  waveguide.  The  VLF  energy  in  the  waveguide 
can  then  leak  out  and  undergo  trans-ionospheric  attenuation,  the  most  important  portion  of  which 
occurs  below  660  km.  In  order  to  properly  validate  the  resulting  power  distribution  with  satellite 
measurements,  propagation  within  the  magnetosphere  must  also  be  considered;  including 
potentially  the  effects  of  the  plasmasphere  and  plasmapause.  A  schematic  for  this  approach  is 
shown  in  Figure  1. 


Figure  1.  Organizational  Framework 

A  number  of  different  levels  of  fidelity  are  possible  for  each  step  in  this  procedure,  and  options 
for  each  of  these  will  be  discussed  in  turn.  Ideally  the  resulting  climatology  would  be  a  function 
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of  frequency,  local  time,  day  of  year,  spatial  location,  as  well  as  some  measure  of  solar  activity, 
and  would  specify  both  the  Poynting  and  wave  normal  vectors.  We  will  attempt  to  meet  these 
requirements  while  being  constrained  by  available  observational  and  computational  resources.  In 
this  initial  study  the  model  will  be  scaled  to  measurements  of  VLF  E-field  spectral  density  taken 
by  the  DEMETER  satellite. 

Some  discussion  of  the  coordinate  systems  employed  here  is  also  warranted.  Wave  data  typically 
organizes  around  magnetic  rather  than  geographic  latitudes,  and  both  systems  will  be  used  here. 
For  magnetic  coordinates  an  apex  system  is  employed  [37],  these  follow  field  lines  and  consist  of 
the  apex  radius  (or  apex  latitude),  apex  longitude,  and  height  above  geoid  (which  specifies  two 
points  on  the  field  line  one  for  the  northern  and  one  for  the  southern  hemisphere).  Field  line 
apexes  were  calculated  using  IGRF1 1  and  a  simple  field  line  tracer  employing  a  fourth  order 
Runga-Kutta  stepping  algorithm  with  a  curvature  dependent  step  size.  One  advantage  of  the  apex 
coordinate  for  our  purposes  is  that  a  given  apex  corresponds  to  a  specific  field  line  and  so  map 
directly  to  the  conjugate  point,  though  there  may  be  a  longitudinal  shift  involved.  Two 
categorizations  of  the  apex  coordinate  will  be  utilized  here:  the  apex  radius  normalized  to  the 
earth's  radius,  A  =  1  +  (hA/Ri),  where  Iia  is  the  maximum  height  (apex)  of  the  field  line  above 
the  geoid,  and  RE  is  the  radius  of  the  earth;  and  the  apex  latitude,  kA  =  ±cos"'(  1/7  A).  Figure  2a 
shows  contours  of  constant  apex  latitude  overlaid  on  a  world  map,  while  Figure  2b  compares  the 
apex  radius  at  a  height  of  700  km  with  the  Mcllwain  L  shell  [38]  at  the  same  height.  Figure  2b 
shows  the  close  correspondence  between  the  Apex  Radius  and  the  Mcllwain  L,  with  average 
contour  values  being  the  same  but  with  the  Mcllwain  L  going  to  higher  values  where  the 
geomagnetic  field  is  relatively  weak  (near  the  south  Atlantic  anomaly)  and  to  lower  values  where 
the  geomagnetic  field  is  relatively  strong,  in  order  to  keep  the  longitudinal  invariant  constant, 

/  =  £*  (1  -  Bt/B)1/2ds\3H\ 


Figure  2.  Coordinate  Systems 


2.2  Wave  Data 

The  wave  data  used  in  this  study  were  obtained  from  the  Instrument  Champ  Electrique 
experiment  [39]  aboard  the  French  micro-satellite  DEMETER.  Originally  launched  in  June  2004 
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at  an  orbit  altitude  of  710  km,  it  was  lowered  to  660  km  in  December  2005.  Here  we  analyze 
electric  field  spectral  amplitudes  in  the  ELF/VLF  range,  from  10  Hz  up  to  20  kHz  at  a  resolution 
of  19.53 125  Hz,  collected  between  January  2005  and  December  2009.  The  stated  sensitivity  and 
dynamic  range  in  the  VLF  are  0.05  pV/m/vHz  and  >80  dB  respectively,  with  an  instrumental 
noise  of  0.03  -  0.05  pV/nWHz.  Only  survey  mode  data  is  used  here  where  the  electric  field 
component  is  perpendicular  to  the  orbital  plane  (sensor  component:  E12).  DEMETER's  orbit  is 
quasi-sun-synchronous  and  so  measurement  times  are  centered  on  two  local  times,  around  1030 
and  2230  for  downward  and  upward  orbits  respectively.  Data  at  latitudes  above  about  65°  were 
not  collected.  We  use  the  “VLF  electric  spectrum”  data  type  (APID  1 132,  Type  0),  which  consist 
of  two  consecutive  spectra  each  individually  consisting  of  an  on-board  average  of  40  spectra 
collected  over  a  duration  of  2.048  seconds.  This  data  product  is  available  nearly  continuously 
during  the  DEMETER  mission. 

In  order  to  investigate  the  data  in  a  climatological  sense  they  are  averaged  into  l°xl0  geographic 
latitude  and  longitude  spatial  bins  for  each  half  orbit  for  each  month  of  the  year.  The  data  are 
also  averaged  over  5  frequency  bins  to  give  a  final  spectral  resolution  of  97.65625  Hz.  The 
averaging  procedure  involves  summing  all  spectral  densities  and  recording  the  number  of 
observations  which  contribute.  The  average  is  then  the  total  spectral  density  divided  by  the  total 
number  of  observations.  This  same  approach  was  taken  when  averaging  across  other  dimensions, 
e.g.  when  averaging  over  apex  latitudes  all  spectral  densities  within  the  range  was  summed  and 
then  divided  by  the  number  of  observations  to  give  an  average.  Although  the  data  is  often 
presented  as  the  log  of  the  spectral  density,  in  no  cases  were  averages  performed  on  the  log 
transformed  data. 


Figure  3.  Histograms  of  Contributing  Observations  from  DEMETER 

Figure  3  shows  a  Histogram  of  the  number  of  observations  contributing  to  each  spatial/temporal 
average  considered  here  during  January  and  July  at  5  kHz.  The  bulk  of  the  locations  have  20  - 
30  observations  making  up  their  averages  (corresponding  to  800  -  1200  individual  spectra),  thus 
giving  us  some  confidence  in  their  validity  as  averages  and  explaining  somewhat  the  binning 
dimensions  employed.  Figure  4  shows  some  representative  spatial  plots  of  the  data  overlaid  on  a 
global  map  of  continental  outlines,  and  Figure  5  shows  the  same  data  plotted  as  histograms.  Both 
Figure  4  and  5  show  averaged  data  from  all  years  of  DEMETER  operation  at  5000  -  5078  Hz 
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separated  into  upward  (nighttime)  and  downward  (daytime)  orbits  for  January  and  July.  Both 


2  2 

figures  show  the  data  in  spectral  density  units  of  log(pVVmVHz)  and  the  histogram  was 
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constructed  using  a  bin  size  of  0.025  log(pV  /m  / Hz). 
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Figure  4.  Representative  VLF  Spectral  Densities  from  DEMETER  Averages 
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2005  -2009,  5.0  -  5.1  kHz 


log(pV2/m2/hlz) 


Figure  5.  Representative  Histograms  of  VLF  Spectral  Densities  from  DEMETER  Averages 

While  the  nighttime  distributions  look  smooth  and  normally  distributed  (in  log(power)  space), 
with  the  mean  value  in  July  being  significantly  higher  than  it  is  in  January,  it  is  clear  that  both 
sets  of  daytime  spatial  averages  push  up  against  the  instrumental  noise  floor  (~  -2.8 
log(pV2/m2/Hz)  at  5  kHz  from  Figure  4  in  Berthelier  et  al.  [39])  in  a  significant  number  of  the 
averaged  values.  It  is  also  noteworthy  that  the  data  does  not  simply  cut  off  at  the  sensitivity 
value,  but  shows  a  second  distribution  as  the  noise  floor  is  approached.  This  reflects  the  fact  that 
each  grid  average  itself  contains  a  distribution  of  events  and  each  spectrum  is  itself  the  on-board 
average  of  40  spectra.  This  has  consequences  for  the  interpretation  of  the  data,  e.g.  the  true  range 
of  magnitudes  is  not  represented  and  unresolved  impulsive  events  are  simply  time  averaged  into 
the  resulting  values.  It  also  accounts  for  the  relatively  wider  peak  associated  with  the  'January 
Day'  values  near  the  noise  floor  as  more  individual  spectra  below  the  cutoff  would  have  been 
integrated  into  each  survey  mode  value. 

In  the  case  of  the  July,  daytime  distribution  it  appears  as  if  there  is  a  symmetric  distribution 
centered  around  -1.5  log(pV  /m“/Hz),  which  is  superimposed  on  the  noise  floor  distribution 
centered  around  -2.8  log(pV2/m2/Hz).  It  is  not  clear  if  a  similar  symmetric  distribution  exists  in 
the  January,  daytime  data,  if  so  it  is  obscured  by  the  noise  floor.  In  fact  it  is  not  immediately 
clear  why  the  entire  spatial  distribution  of  log(power)  VLF  should  be  normally  distributed  at  all. 
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This  implies  that  the  distribution  is  created  by  a  series  of  uncorrelated  events  whose  effects  are 
multiplicative;  perhaps  the  distribution  of  power  from  lightning  strokes  and  the  absorption  due  to 
transionospheric  propagation.  Their  covariance  should  be  small  but  their  contribution  to  power 
measured  at  the  satellite  is  multiplicative  (Sferic  Power  x  Transionospheric  Absorption  oc  Power 
Received). 

It  is  of  note  that  both  Figure  4  and  5  show  significant  diurnal  and  seasonal  variation.  These 
features  along  with  systematic  spatial  features  will  be  discussed  in  later  sections.  As  a  check  on 
the  year  to  year  uniformity  of  the  data  an  additional  histogram  of  the  July  averages  at  5  kHz  for 
each  year  of  data  individually  is  presented  in  Figure  6.  The  shape  of  the  distribution  function  is 
nearly  identical  within  each  year  with  variation  mostly  in  the  total  number  of  samples  collected, 
e.g.  there  were  significantly  fewer  samples  archived  in  2005,  the  first  full  year  of  the  DEMETER 
mission.  There  is  also  a  slight  preponderance  of  values  above  1  pV  / m  /Hz  in  the  daytime  2005 
values  compared  with  values  between  0.01  and  0. 1  pV2/m2/Hz,  this  is  most  likely  due  to  the 
more  active  sun  during  that  period.  These  emissions  belong  to  class  of  intense  typically 
broadband  emission  observed  along  the  auroral  oval  (as  seen  in  Figure  4)  which  don't  show  the 
seasonality  of  lightning.  It  is  not  expected  that  the  change  in  altitude  from  -700  to  -660  km  at 
the  end  of  2005  would  have  a  major  impact  on  transionospheric  attenuation  since  the  majority  of 
absorption  takes  place  at  the  lower  edge  of  the  ionosphere  (60-100  km)  where  ion-neutral 
collisions  are  high,  and  the  remaining  part  is  due  to  Coulomb  interactions  integrated  along  the 
path  [40]  which  should  peak  at  the  F  layer.  It  is  likely  that  these  emissions  are  related  to  events  in 
the  outer  magnetosphere  which  are  more  easily  driven  by  solar  conditions.  It  is  also  possible  that 
more  active  solar  conditions  would  affect  amplification  of  whistler  energy  during  its  passage 
through  the  inner  magnetosphere.  In  our  analysis  of  lightning  we  will  in  general  take  measures  to 
exclude  auroral  phenomena. 
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July,  5.0  -  5.1  kHz 


log((.iV2/m2/Hz) 

Figure  6.  Representative  Annual  Histograms  of  VLF  Spectral  Densities  from  DEMETER 
2.3  Optical  Lightning  Data 

The  optical  lightning  flash  data  used  here  were  collected  by  the  Optical  Transient  Detector 
(OTD)  and  its  successor  instrument  the  Lightning  Imaging  Sensor  (LIS).  Inter-calibrated  and 
merged  data  products  were  downloaded  from  a  NASA  maintained  website 
(http://thunder.msfc.nasa.gov/data/index.html)  between  October  and  December  2010.  Two 
different  available  climatologies  were  combined  for  our  purposes:  the  HRMCCOMFR 
(monthly),  which  has  a  0.5  degree  spatial  resolution  and  monthly  temporal  resolution;  and  the 
LRADCCOMSMFR2,  which  has  2.5  degree  spatial  resolution  and  2  hr  temporal  resolution  for 
every  day  of  the  year.  From  these  a  combined  climatology  with  1  degree  spatial  resolution  and 
an  average  diurnal  cycle  (with  2  hr  resolution)  for  each  average  monthly  value  were  constructed. 
The  averaging  procedure  involved  deresolving  the  HRMC  COM  FR  to  1  degree  spatial 
resolution,  then  constructing  an  average  diurnal  cycle  for  each  month  from 
LRADC  COM  SMFR2  and  interpolating  that  to  1  degree  resolution,  this  was  normalized  and 
applied  to  the  deresolved  HRMC  COM  FR  set.  Figure  7a  shows  the  global  monthly  averages 
and  Figure  7b  the  global  annual  average  diurnal  cycle  of  the  combined  climatology.  The  two 
vertical  lines  in  Figure  7b  represent  the  average  local  time  of  DEMETER  overpasses.  Both 
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Figure  7a  and  7b  show  the  canonical  average  of  44±5  lightning  flashes  occurring  around  the 
globe  every  second  [41]. 


Figure  7.  Temporal  Optical  Flash  Climatology 

Figure  7a  also  shows  the  monthly  averages  separated  into  hemispheres  (red  for  the  Southern 
Hemisphere  and  blue  for  the  Northern  Hemisphere),  and  the  dominance  of  the  northern 
hemisphere  on  the  overall  global  cycle  is  clearly  evident,  this  effect  is  primarily  due  to  the 
preponderance  of  land  mass  in  the  northern  hemisphere.  Figure  7b  also  shows  the  diurnal  cycle 
broken  into  the  major  land  mass  areas  (along  with  surrounding  oceans  following  the  grouping  of 
Lay  et  ah,  [42])  as  different  colored  lines  (North  America  in  blue,  South  America  in  light  blue, 
Europe  in  green,  Africa  in  light  green,  Asia  in  yellow,  and  Australia  in  red)  and  while  we  see  a 
universality  to  the  shape  of  the  diurnal  cycle  around  the  world  there  are  also  variations  in  peak 
time  and  relative  amplitude.  The  region  to  region  differences  shown  here  are  consistent  with 
those  shown  at  a  higher  temporal  resolution  in  Mach  et  al.  [43]  (also  from  OTD/LIS  data).  Mach 
et  al.  [43]  also  present  OTD/LIS  data  limited  to  oceanic  observations.  These  show  a  much  less 
pronounced  diurnal  cycle  with  a  maximum  in  the  early  morning  when  land  based  lightning  is 
near  its  minimum.  These  features  are  also  extant  in  the  combined  climatology  constructed  here. 

Lay  et  al.  [42]  presented  diurnal  lightning  statistics  from  World  Wide  Lightning  Location 
Network  (WWLLN)  data,  and  these  showed  a  larger  range  of  regional  peak  times  and 
amplitudes.  A  recent  comparison  between  WWLLN  and  the  National  Lightning  Detection 
Network  (NLDN)  over  North  America  suggests  that  this  difference  stems  partially  from 
WWLLN’s  bias  towards  stronger  flashes,  partially  from  NLDN’s  limited  sampling  of  lightning 
over  the  ocean,  and  remains  partially  other/unexplained  [44].  In  regions  of  overlap  the  NLDN 
diurnal  cycle  agrees  closely  with  the  OTD/LIS  data.  The  global  diurnal  peak  in  lightning  occurs 
in  the  late  afternoon  in  all  three  data  sets.  The  local  time  associated  with  DEMETER  daytime 
overpasses  is  near  the  annually  averaged  diurnal  minimum  such  that  more  lightning  would  be 
expected  during  nighttime  DEMETER  observations.  This  is  contrary  to  tautology  which  puts 
more  lightning  during  the  day.  Figure  8  shows  some  representative  spatial  plots  of  the  combined 
climatology,  for  January  and  July  200  -  400  and  1400  -  1600  hours  local  time.  The  minimum 
flash  rate  shown  has  been  artificially  limited  at  10'4  times  the  maximum  flash  rate  to  enhance 
details.  The  preponderance  of  lightning  flashes  over  the  land  is  immediately  evident  with  the 
average  land  to  ocean  ratio  being  about  a  factor  of  ten.  The  effects  of  season  is  also  evident  with 
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the  summer  hemisphere  showing  lightning  extending  to  higher  latitudes.  The  features  evident  in 
Figure  8  are  consistent  with  previous  publications  of  the  OTD/LIS  data,  e.g.  [41]. 


Flashes/km2/day,  January,  LT  =  14  -16 


Flashes/km2/day,  July,  LT  =  14  -16 


80 

60 

40 

20 

0 

-20 

-40 

-60 

-80 


-160  -120  -80  -40  0  40  80  120  160 


80 

60 

40 

20 

0 

-20 

-40 

-60 

-80 


-160  - 


120 


-80 


-40 


120 


160 


c 


Flashes/km2/day,  January,  LT  =  2-4 


80 

60 

40 

20 

0 

-20 

-40 

-60 

-80 


-160  -120  -80  -40  0  40  80  120  160 


-160  -120  -80  -40  0  40  80  120  160 


Flashes/km2/day,  July,  LT  = 


40 

20 

0 

-20 

-40 

-60 

-80 


Figure  8.  Representative  Spatial  Optical  Flash  Climatology 


2.4  Interpreting  the  Wave  Data 

In  order  to  constrain  our  initial  model  with  measurements  of  VLF  spectral  density  taken  by  the 
DEMETER  satellite,  the  contribution  of  lightning  to  the  DEMETER  data  set  must  be  assessed 
and  organized.  The  temporal  and  spatial  dependencies  of  VLF  from  LGW  in  LEO  are  expected 
to  be  different  from  those  of  space  based  sources  of  VLF  radiation  such  as  plasmaspheric  hiss. 
Lightning  occurs  predominantly  over  land,  in  the  evening,  and  during  summer.  These  are 
characteristics  associated  with  terrestrial  meteorology  and  are  not  expected  to  be  shared  by  non¬ 
terrestrial  VLF  sources  observed  in  LEO.  In  addition  to  considering  the  climatology  of  lightning 
in  interpreting  DEMETER  measurements  it  must  be  considered  that  the  VLF  from  lightning 
propagates  up  through  the  ionosphere  before  being  detected.  Conversion  of  VLF  radiation  in  the 
earth  ionosphere  waveguide  into  whistler  mode  radiation  beyond  the  ionosphere  is  a  complicated 
process  and  depends  in  a  non- trivial  way  on:  the  total  electron  and  ion  content  and  profile  in  the 
ionosphere,  the  dip  angle  and  intensity  of  the  geomagnetic  field,  the  frequency  of  the  incident 
radiation,  and  even  the  wave-normal  of  the  radiation  with  respect  to  the  ionospheric  electron 
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gradient  and  geomagnetic  field.  Still,  as  a  general  rule  there  is  more  absorption  at  low  magnetic 
latitudes,  during  the  day,  and  at  higher  frequencies.  These  relationships  are  typically 
characterized  by  power  laws  [40].  It  can  be  clearly  seen  in  Figure  4a  that  the  minimum  VLF 
spectral  densities  at  DEMETER  altitudes  occur  right  along  the  magnetic  equator  and  not  the 
solar  equator  as  observed  for  lightning.  We  know  that  VLF  trajectories  in  magnetized  plasma 
follow  field  lines  directly  or  indirectly.  These  observations  confirm  the  importance  of  LGW 
conversion  and  attenuation  in  understanding  the  DEMETER  observations  and  make  the  use  of 
both  magnetic  and  geographic  coordinates  essential. 

Globally  lightning  peaks  in  July,  so  we  expect  VLF  power  from  lightning  to  also  show  a  peak  in 
July.  Looking  for  hemispheric  trends,  as  indicated  by  Figure  7,  is  more  complicated  because 
ducted  whistler  radiation  also  shows  up  at  the  conjugate  location,  and  even  in  the  absence  of  a 
duct  whistler  propagation  is  strongly  affected  by  the  geomagnetic  field  and  plasmapause  location 
such  that  power  may  show  up  in  the  conjugate  hemisphere  at  a  similar  apex  latitude.  As  a  first 
cut  Figure  9a  shows  a  plot  of  pV7m  /Hz  vs.  frequency  for  nighttime  orbits  and  aa  below  55°, 
Figure  9b  shows  daytime  orbits  and  aa  below  48°.  Both  figures  show  the  data  with  monthly 
averages  plotted  in  different  colors  where  violet  -t>  red  correspond  to  January  -t>  December. 
Higher  Apex  latitudes  have  been  excluded  because  of  strong  non-lightning  emissions  there 
associated  with  auroral  processes;  also,  as  seen  in  Figure  8,  the  majority  of  lightning  occurs  at 
lower  latitudes.  The  reason  for  the  relatively  lower  latitude  cutoff  during  the  day  is  due  to  a 
combination  of,  diurnal  variation  of  the  auroral  oval  latitude  and  the  relatively  weaker  spectral 
amplitude  during  the  DEMETER  daytime  orbits  making  any  contamination  more  problematic.  In 
Figure  9a  the  blue  trace  representing  March  still  shows  features  around  6-8  kHz  possibly 
associated  with  high  latitude  processes. 


0  5  10  15  20  0  5  10  15  20 

Frequency  (kHz)  Frequency  (kHz) 

Figure  9.  DEMETER  Spectral  Density  vs.  Frequency 

There  is  a  clear  annual  cycle  evident  in  Figure  9  with  maximum  power  being  seen  in  July  and 
August  as  expected  for  a  lightning  driven  process.  This  annual  variation  is  seen  in  both  daytime 
and  nighttime  orbits  with  similar  amplitude.  The  nighttime  measurements  show  higher  values  by 
an  order  of  magnitude  or  more,  as  expected  given  the  larger  occurrence  of  lightning,  the  higher 
altitude  bottom  of  the  ionosphere,  and  the  smaller  total  electron  content  expected  during 
DEMETER  upward  orbits.  Both  plots  show  a  broad  peak  around  3-5  kHz  consistent  with 
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lightning  after  undergoing  frequency  dependent  absorption  in  the  ionosphere.  The  sharp  peaks  in 
the  spectra  are  due  to  VLF  transmitters  and  they  show  similar  diurnal  variation,  but  no  additional 
annual  cycle.  The  daytime  orbits  show  a  maximum  below  1  kHz  with  little  annual  cycle  which  is 
consistent  with  a  Hiss  like  emission.  The  nighttime  feature  in  this  frequency  range  does  show  the 
same  annual  cycle  and  so  is  likely  driven  by  lightning.  The  flatness  of  the  spectrum  is  surprising 
especially  for  the  nighttime  data  and  it  is  clear  that  the  seasonality  extends  to  high  frequencies.  If 
binned  more  finely  these  spectral  features  do  show  systematic  latitudinal  trends,  but  do  not  seem 
to  depend  strongly  on  proximity  to  an  area  of  active  lightning.  Some  additional  discussion  of 
these  features  is  included  below  and  in  later  sections  but  a  full  analysis  is  beyond  the  scope  of  the 
current  investigation. 


A  cursory  inspection  of  figures  4  and  8  show  the  general  spatial  correspondence  of  the 
DEMETER  VLF  data  at  5  kHz  with  the  OTD/LIS  lightning  maps  below  auroral  latitudes,  and 
Figure  9  confirms  that  the  temporal  variation  is  consistent  with  lightning  as  well.  In  Figure  10  a 
color  plot  of  VLF  spectral  density  is  presented,  again  excluding  high  latitudes,  as  a  function  of 
longitude  and  frequency.  The  purpose  here  is  to  break  down  the  range  of  frequencies  which  show 
a  spectral  and  spatial  relationship  consistent  with  lightning.  The  months  of  January  and  July  are 
presented  separated  into  upward  and  downward  orbits. 
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Figure  10.  DEMETER  Spectral  Density,  Frequency  vs.  Longitude 


14 


In  Figure  10  a-d  VLF  transmitters  show  up  as  horizontal  bands,  with  NWC  being  most 
prominent  at  a  frequency  of  19.8  kHz  and  1 14°  longitude  (it  also  shows  the  spectral  broadening 
observed  previously,  e.g.  Parrot  et  ah,  [45]).  The  Hiss  like  emission  below  1  kHz  seen  in  Figure 
9  shows  up  at  all  longitudes  during  the  day.  Overall  the  most  evident  features  of  Figure  10  show 
qualitative  agreement  with  the  spatial  and  temporal  variation  associated  with  lightning.  More 
VLF  spectral  density  is  present  during  nighttime  orbits,  and  in  the  northern  hemisphere  summer. 
There  are  major  features  associated  with  land  masses  in  their  summer  hemispheres  consistent 
with  Figure  3.  Figure  10a  (January  daytime  orbits)  shows  the  lowest  magnitudes  and  very  little  in 
the  way  of  systematic  patterns,  with  the  exception  of  the  Hiss  band  already  mentioned.  Referring 
back  to  Figure  4a  we  see  that  a  large  fraction  of  the  observations  in  Figure  10a  are  showing  the 
effects  of  instrumental  sensitivity  and  these  have  been  averaged  into  the  results  shown.  Figure  10 
along  with  Figure  4  make  it  clear  that  the  broad  peak  at  5  kHz  seen  in  Figures  9a  and  9b  are 
strongly  correlated  to  major  landmasses  and  so  lightning.  This  is  also  consistent  with  the 
seasonal  variation  shown  there  and  the  spectral  content  of  lightning.  The  nighttime  peak  near  15 
kHz  shown  in  Figure  9b  also  correlates  with  the  landmasses  and  shows  an  interesting  frequency 
dependence  as  a  function  of  longitude. 


In  order  to  investigate  these  features  in  more  detail  Figures  1 1  and  12  show  similar  plots 
(spectral  densities  for  July  and  January  respectively  separated  into  upward  and  downward  orbits 
as  a  function  of  frequency  and  longitude)  now  segregated  according  to  magnetic  hemisphere 
(apex  latitudes  between  15  and  50  for  the  Northern  Hemisphere  and  -15  to  -50  for  the  Southern 
Hemisphere).  In  these  figures  tropical  latitudes  have  been  mostly  eliminated,  and  the  scale  has 
been  adjusted  in  order  to  emphasize  seasonal  trends.  Also  included  in  Figures  1 1  and  12  are  the 
averaged  lower  hybrid  resonance  frequencies  (LHR)  calculated  from  averaged  plasma 
compositions  estimated  with  the  International  Reference  Ionosphere  (IRI)  combined  with 
averaged  geomagnetic  field  values  from  the  International  Geomagnetic  Reference  Field  (IGRF). 
The  most  updated  versions  of  these  reference  codes  at  the  time  of  calculation  were  utilized  (IRI- 
2007  and  IGRF-201 1).  The  local  ion  composition  for  2007  was  calculated  with  the  IRI  and  the 
local  geomagnetic  field  for  2007  was  calculated  with  the  IGRF  (2007  was  chosen  because  it  was 
the  latest  year  available  for  the  IRI  and  it  was  the  average  year  of  the  DEMETER  data  utilized), 
the  middle  day  of  the  month  was  chosen  for  the  local  times  relevant  to  the  DEMETER 

overpasses.  The  calculation  of  the  LHR  frequency  strictly  follows  that  described  in  Brice  and 

2  2 
1  LJr  LJlf 

Smith  [46],  =  ^  ,  a  . ,  "5  where  Meff  is  the  plasma  effective  mass,  cop  is  the  electron 

Msff 


plasma  frequency,  and  coh  is  the  electron  cyclotron  frequency. 
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Figure  11.  Hemispheric  July  Spectral  Densities,  Frequency  vs.  Longitude 
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Figure  12.  Hemispheric  January  Spectral  Densities,  Frequency  vs.  Longitude 


Figures  1 1  and  12  again  show  major  features  consistent  with  lightning  and  present  an 
explanation  for  the  higher  frequency  waves  showing  an  annual  cycle  in  Figures  9a  and  9b.  They 
are  consistent  with  the  LHR  Hiss  described  by  Brice  and  Smith  [46]  and  observed  innumerable 
times  since.  In  that  work  they  are  described  as  being  generated  in  the  immediate  vicinity  of  the 
satellite,  associated  with  nonducted  whistler  propagation  from  the  opposite  hemisphere,  and 
exhibiting  a  lower  cutoff  at  the  LHR.  This  lower  cutoff  reflects  the  fact  that  the  LHR  its  self 
defines  a  cutoff  frequency  for  propagation  transverse  to  the  ambient  magnetic  field.  Such  quasi¬ 
resonant,  quasi-electrostatic,  whistlers  can  show  a  large  increase  in  intensity  close  to  the  LHR  as 
their  group  velocity  decreases  and  they  near  magnetic  reflection.  Additional  discussion  of 
relevant  LHR  observations  and  theoretical  treatments  can  be  found  in  e.g.  [21,  47,  48,  49]. 
Figures  1 1  and  12  show  that  the  LHR  emissions  are  associated  with  lightning  in  the  conjugate 
hemisphere,  so  that  Figure  1 1  shows  maxima  above  the  LHR  in  the  southern  hemisphere 
associated  with  lightning  over  North  America  and  Euro/ Asia  in  the  northern  hemispheric 
summer.  The  shift  in  the  predicted  lower  cutoff  of  the  LHR  hiss  band  between  day  and  night  is 
consistent  with  the  calculated  values.  Binning  the  results  into  smaller  latitudinal  regions  results 
in  good  agreement  between  the  shifting  LHR  and  the  lower  cutoff  frequency  (not  shown).  In  this 
first  cut,  single  wave-normal  component,  climatological  look  at  the  data  it  is  not  possible  to 
distinguish  uniquely  between  the  local  quasi-resonant  LHR  hiss  and  a  propagating  quasi- 
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longitudinal  wave  above  the  LHR  (e.g.  a  0+  whistler,  or  a  ducted  1-  whistler).  Knowledge  of  the 
total  wave  nonnal  vector  would  potentially  allow  discrimination,  or  alternatively  reliable 
observations  of  the  magnetic  field  as  the  quasi-electrostatic  waves  have  very  little  power  in  their 
AC  magnetic  fields. 

Figures  9-12  show  that  the  major  features  of  the  DEMETER  VLF  data  set  above  1-2  kHz  can  be 
explained  either  by  lightning  or  auroral  processes.  This  conclusion  is  consistent  with  previous 
work  by  Nemec  et  al.  [50].  The  spatial,  annual,  and  diurnal  features  have  all  been  explored  and 
along  with  spectral  properties  are  found  to  be  consistent  with  lightning.  Although  only  one  wave 
nonnal  component  is  represented  in  the  data  theoretical  predictions  and  satellite  observations 
confirm  that  the  wave  vector  is  near  vertical  within  the  ionosphere  e.g.  [51],  at  least  for  0+ 
whistlers,  and  so  we  expect  the  measured  component  to  capture  most  of  the  field  intensity.  In  any 
case  in  concert  with  the  IGRF  or  the  satellite  ephemeris  we  know  the  angle  between  the  wave 
vector  and  the  local  geomagnetic  field  and  so  could  reconstruct  the  full  intensities  from  the 
polarization  ratios.  Nemec  et  al.  [50]  also  concluded  that  DEMETER  survey  mode  data 
“represents  a  sufficient  proxy  for  the  intensity  of  electromagnetic  waves”.  This  data  set  then, 
provides  an  excellent  benchmark  for  a  climatology  of  lightning  generated  VLF  in  LEO 
constructed  by  other  means. 

2.5  VLF  from  Lightning  in  the  Earth-Ionosphere  Waveguide 

The  lightning  climatology  employed  here  has  already  been  outlined  in  section  4,  it  originates 
from  impulsive  optical  measurements  taken  from  space  and  can  be  specified  in  flashes/km2/day 
as  shown  in  Figure  7.  Our  objective  is  to  construct  a  climatology  of  VLF  radiation  from  lightning 
in  low  earth  orbit,  so  we  would  like  to  convert  flashes/km2/day  in  the  troposphere  to 
Power/km  /Hz  at  LEO.  The  initial  approach  taken  here  assumes  that  the  number  of  flashes/km" 
is  proportional  to  the  amount  of  power/km2  radiated  in  the  VLF.  This  is  equivalent  to  saying  that 
the  strikes  are  incoherent  and  the  distribution  of  power  over  all  strikes  is  the  same  as  a  function 
of  time  and  place.  There  is  some  evidence  to  suggest  that  this  is  not  quite  the  case  as,  for 
example,  the  diurnal  cycle  of  flash  number  is  quite  different  over  oceans  [43]  and  so  might  be  the 
power  distribution  [52].  Initial  studies  of  seasonal  variation  in  the  charge  moment  distribution 
functions  for  intense  lightning  indicate  overall  similarity  to  the  total  [35].  However,  there  are 
also  suggestions  that  the  temporal/spatial  variation  of  powerful  strokes  is  somewhat  different 
than  that  of  the  entire  range  of  powers  as  touched  on  in  Section  4  and  alluded  to  in  Lay  et  al. 

[42].  While  it  would  certainly  be  of  interest  to  include  a  spatially  and  temporally  dependent  flash 
to  power  conversion  factor  there  is  insufficient  guidance  on  this  point  in  the  literature  to  do  so  at 
this  time. 

Once  the  spatial  and  temporal  distribution  of  VLF  power  generation  from  lightning  has  been 
fixed  it  is  still  necessary  to  propagate  that  power  within  the  earth-ionosphere  waveguide  and  then 
couple  it  into  whistler  mode  radiation  at  the  bottom  of  the  ionosphere.  At  the  most  general  level 
we  expect  power  spreading  in  free  space  to  drop  off  as  the  square  of  the  distance  due  to  purely 
geometric  considerations.  Because  the  conductive  ground  and  conductive  ionosphere  form  a 
spherical  waveguide  for  VLF  radiation,  at  the  next  level  of  abstraction  we  expect  that  attenuation 
will  transition  to  an  inverse  distance  proportionality  as  distance  increases,  and  show  a 
pronounced  frequency  dependence.  At  even  longer  distances  we  expect  to  see  some  focusing 
occur  because  the  earth  is  round  and  rays  diverging  from  an  initial  point  will  start  to  converge  at 


18 


the  antipodal  point.  The  exact  details  of  how  power  in  the  waveguide  couples  to  the  whistler 
mode  are  still  obscure,  and  preliminary  studies  suggest  that  modeled  electric  fields  are  more  than 
10  times  smaller  than  in  situ  measurements  in  the  lower  ionosphere  [53], 

Rocket  and  satellite  observations  show  that  sferics  can  couple  to  the  whistler  mode  over 
distances  of  at  least  1000  km  from  the  flash  [32,  33,  51].  In  the  most  complete  study  of  its  kind 
Fiser  et  al.  [54]  analyzed  -30000  lightning  whistler  pairs  and  found  that  the  average  amplitude  of 
a  whistler  is  largest  when  the  magnetic  footprint  of  the  detecting  satellite  is  approximately  1 
degree  north  of  the  causative  lightning  location.  The  spatial  resolution  of  this  study  was  1  degree 
in  latitude  and  all  comparisons  were  made  in  Europe.  This  distance  is  significantly  different  than 
the  vertical  projection  of  the  satellite  location,  which  would  be  about  4  degrees  north  of  the 
magnetic  projection  at  the  location  in  question.  For  distances  in  the  range  of  300-2000  km  they 
found  that  the  mean  whistler  amplitude  in  LEO  drops  off  in  inverse  proportion  to  the  distance 
from  the  location  of  maximum  amplitude.  In  light  of  these  studies  we  utilize  an  inverse  square 
relationship  between  a  flash  location  and  the  power  it  delivers  to  the  ionosphere.  A  more 
complete  treatment  might  transition  from  the  inverse  square  treatment  utilized  here  to  the 
dB/Mm  treatment  more  familiar  to  radio  engineers  at  some  distance  from  the  flash. 

In  the  future  it  would  be  of  interest  to  study  the  frequency  dependence  of  whistler  energy  as  a 
function  of  distance  from  the  flash  location.  Ferencz  et  al.  [55]  report  the  observation  of  “Spiky 
Whistlers”  in  LEO  which  they  attribute  to  the  formation  of  a  whistler  from  a  tweek,  a  tweek 
being  a  sferic  which  has  traveled  far  enough  in  the  earth  ionosphere  waveguide  to  become 
dispersed  at  a  frequency  consistent  with  the  cutoff  frequency  of  the  waveguide,  1.5-2  kHz,  and 
its  hannonics.  This  demonstrates  that  the  frequency  characteristics  of  the  generating  sferic  can  be 
transferred  onto  the  whistler.  It  is  well  known  that  the  frequency  characteristics  of  sferics  change 
greatly  as  a  function  of  distance  from  the  flash.  Such  changes  are  often  described  in  the  context 
of  a  mode  theory  for  waveguide  transmission,  as  described  in  [56].  To  summarize:  frequencies 
around  and  below  5  kHz  are  attenuated  significantly  relative  to  higher  frequencies  which  then 
appear  enhanced  in  the  resulting  spectrum.  During  the  night  this  relative  maximum  in  attenuation 
shifts  somewhat  lower  in  frequency  due  to  the  taller  waveguide,  and  overall  attenuation  is  lower 
due  to  the  smaller  ion/neutral  collision  frequency  at  that  higher  altitude.  Sferics  which  cross  the 
day/night  terminator  will  experience  an  additional  loss  due  to  the  abrupt  change  in  waveguide 
height.  Measurements  taken  nearer  to  a  flash,  and/or  specifically  designed  to  measure  it's 
incipient  spectrum,  show  a  peak  around  3-5  kHz  and  a  drop  off  well  described  by  one  over  the 
frequency  out  into  the  MHz  range  [57].  These  two  competing  characteristic  (incipient  spectrum 
and  transport)  lead  to  a  typical  measured  spectrum  in  the  earth-ionosphere  waveguide  which 
peaks  around  10-15  kHz  (Dowden,  Brundell,  and  Rodger  [58]  suggest  that  half  power 
frequencies  of  5  and  19  kHz  are  reasonable).  Spectra  of  this  type  are  readily  available  on  the 
web,  including  many  associated  with  the  WWLLN  network  (webflash.ess.washington.edu).  This 
description  is  also  borne  out  by  a  recent  publication  on  the  development  of  a  typical  wavefonn 
bank  for  use  in  the  VLF  geolocation  of  lightning  [34]  which  incorporates  thousands  of  measured 
spectra  at  various  distances.  The  Fourier  transform  of  Figure  3  in  that  paper  (which  shows 
expected  received  waveforms  as  a  function  of  distance  from  the  stroke)  would  give  typical 
spectrum  as  a  function  of  distance  from  the  flash  consistent  with  the  above  expectations  (Said, 
private  correspondence). 
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2  2 

In  this  work  the  conversion  from  flashes/km  /day  to  power/km  /day  proceeds  as  follows.  For  a 
given  location  it  is  first  assumed  that  the  number  of  flashes/km2/day  converts  directly  into  a 
power/km  / day  at  that  location  (Po  in  equation  1).  Then  contributions  from  every  other  location 
are  calculated  based  on  a  reference  distance  reflecting  the  length  scale  of  the  contributing  grid 
location  (geocentric  coordinates  are  not  equal  area)  and  a  great  circle  distance  between  the  grid 
centers,  an  additional  normalization  is  included  to  account  for  the  difference  in  areas  between  the 
two  locations  and  so  conserve  the  per  km~  designation. 


ref  erence  distance *  area(i) 

-  * -  (1) 

greatcdistance  area(O) 

t=i 

A  limit  on  the  great  circle  distance  which  can  contribute  to  a  location  has  not  been  included  at 
this  time  nor  has  antipodal  focusing.  This  calculation  is  done  in  universal  time  and  then 
converted  to  local  time.  Figure  13  shows  representative  plots  of  the  output  which  we  refer  to  as 
pseudopower  for  January  and  July  1 100  and  2300  hours  local  time.  The  units  are  arbitrary  but 
the  calculations  were  done  conserving  relative  values  per  day  and  per  km  .  In  Figure  13  we  have 
multiplied  by  the  grid  area  at  each  location  (giving  units  of  pseudopower/day)  in  order  to 
emphasize  features  related  to  the  spatial  distribution  of  power  rather  than  the  unequal  area  of  our 
geographic  grid  which  gets  quite  small  at  high  latitudes. 


20 


Figure  13.  Representative  Spatial  Pseudopower  Climatology 


Attenuation  of  VLF  energy  in  the  earth  ionosphere  waveguide  is  often  discussed  in  terms  of 
dB/1000  km,  e.g.  [56,  59]  with  typical  values  being  a  few  dB/Mm.  Our  inverse  square  treatment 
with  reference  distances  around  100  km  (1°  geographic)  drops  off  very  rapidly  at  first,  a  factor  of 
100  in  the  first  1000  km  corresponding  to  20  dB/Mm,  then  more  slowly  equivalent  to  2.5  dB/Mm 
in  the  4th  1000  km.  Over  the  first  10  Mm  the  intensity  drops  off  by  10"4  corresponding  to  an 
average  attenuation  of  4  dB/Mm.  At  longer  distances  the  attenuation  becomes  even  smaller  ~1 
db/Mm  in  the  10th  1000  km.  So,  while  simplistic,  the  inverse  square  treatment  provides  a  good 
first  cut  and  for  our  purposes  is  superior  to  a  pure  dB/Mm  approach. 

As  noted  above  the  initial  spectrum  of  lightning  is  thought  to  peak  around  3-5  kHz  and  then 
drop  off  as  one  over  the  frequency  out  to  at  least  the  20  kHz  considered  here.  Considering  that 
the  peak  spectrum  moves  to  higher  frequencies  with  distance,  and  that  empirically  the  power  of  a 
satellite  detected  0+  whistler  drops  off  as  I  /distance”,  it  is  reasonable  to  adopt  a  generic  spectrum 
shape  extracted  from  Lauben,  Inan,  and  Bell  [24]  and  given  in  Equation  2,  which  shows  a  broad 
peak  between  2  and  6  kHz.  In  Equation  2,  to  is  the  angular  frequency  while  a  and  b  are  constants 
(here  a  =  5000  and  b  =  1000).  It  is  normalized  over  the  frequency  range  of  interest  (0  -  20  kHz) 
before  application  to  the  pseudopower  distributions. 
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2.6  Transionospheric  Attenuation 

Once  we  have  a  VLF  power  spectrum  in  the  earth  ionosphere  waveguide  we  must  approximate 
in  some  way  the  leakage  into  the  ionosphere  and  the  attenuation  as  the  wave  energy  propagates 
in  the  plasma.  Here  we  will  use  the  Helliwell  curves  (figure  3-35  in  [40]),  which  are  calculated 
employing  refractive  indices  from  the  quasi -longitudinal  approximation  to  Appleton's  equation. 
Helliwell's  calculations  assumed  a  specific  form  for  the  average  electron  density  and  collision 
frequency  profiles  in  the  ionosphere,  for  daytime  and  nighttime  conditions  separately.  Then  the 
path  integrated  absorption  was  calculated  through  the  ionosphere  for  vertically  incident  radiation 
of  2  and  20  kHz  at  a  range  of  geomagnetic  latitudes  above  10°.  In  our  application  absorption  for 
intennediate  frequencies  are  estimated  assuming  that  absorption  is  proportional  to  the  square  root 
of  the  frequency. 

The  validity  of  the  Helliwell  approach  is  checked  in  some  detail  by  Tao,  Bortnik,  and  Friedrich 
[60],  who  find  reasonable  agreement  with  a  more  numerically  intensive  full  wave  calculation 
coupled  to  more  modem  values  for  the  electron  density  and  collision  frequency  profiles. 

Lehtinen  and  Inan  [61]  also  present  state  of  the  art  calculations  which  are  compared  to  other 
Helliwell  based  calculations  perfonned  by  other  researchers.  The  main  limitations  associated 
with  application  of  the  Helliwell  curves  are  twofold:  first,  the  quasi-longitudinal  approximation 
upon  which  they  are  based  is  known  to  fail  for  large  values  of  the  geomagnetic  inclination  angle 
such  as  exist  near  the  equator;  and  second,  the  ionospheric  profiles  employed,  while  reasonably 
representative,  do  not  show  the  range  of  systematic  variation  now  known  to  exist  in  time  and 
space  (as  summarized  by  e.g.  IRI).  Nevertheless  the  Helliwell  curves  provide  a  widely  accepted 
baseline  from  which  to  proceed.  The  difficulty  of  using  updated  absorption  values  is  due  in  part 
to  the  wide  variety  of  variables  that  transionospheric  absorption  will  depend  on,  including 
geomagnetic  inclination  angle,  geomagnetic  field  strength,  electron  density  and  collision  profile. 
Each  of  which  have  their  own  spatial  and  temporal  dependencies.  Thus  in  our  application  we 
would  like  to  have  an  average  absorption  calculated  for  climatological  averages  of  each  two 
hours  in  local  time,  each  month,  each  frequency,  at  each  spatial  location. 
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2  kHz,  Winter,  1200  hours  Universal  Time 

10  100 
Absorption  (dB) 


2  kHz,  Summer,  0  hours  Universal  Time 


Figure  14.  Transionospheric  Absorption,  Representative  Spatial  Climatology 

In  our  particular  application  the  absolute  values  of  the  absorption  are  not  as  important  as  the 
relative  values,  and  so  details  such  as  polarization  loss  (3  dB),  reflection  (1-2  dB),  and  whistler 
coupling  factor  (discussed  in  e.g.  [24])  while  potentially  important  are  not  examined  in  any 
detail.  Our  values  are  taken  directly  from  Helliwell  [40]  Figure  3-35.  Figure  14  shows  some 
representative  global  transionospheric  absorption  plots  at  2  and  20  kHz  for  winter  and  summer  at 
noon  and  midnight.  Temporal  transitions  are  handled  heuristically  in  a  manner  similar  to  Starks 
et  al.  [18],  involving  an  inverse  tangent  function  of  the  solar  zenith  angle  at  D  region  altitudes.  In 
local  time  coordinates  there  is  no  longitudinal  dependence  to  the  Helliwell  curves  and  so,  for 
example,  the  appropriate  curve  to  apply  to  Figure  13a  at  2  kHz  would  be  a  longitudinal  section  of 
Figure  14a  at  -15°  longitude,  and  13c  at  20  kHz  would  take  a  longitudinal  section  of  figure  14c  at 
165°  longitude. 
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PseudoPower/km2/doy,  January,  1100  hours,  5.0  -  5.1  kHz 
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Figure  15.  Pseudopower  in  Low  Earth  Orbit,  Representative  Spatial  Climatology 


Figure  15  shows  the  results  of  applying  the  appropriate  Helliwell  absorptions  at  5  kHz  to  the 
pseudopower  distributions  shown  in  Figure  13.  The  values  have  been  normalized  to  a  maximum 
of  100  and  plotted  as  log(pseudopower/km2/day)  to  facilitate  comparison  with  Figure  3. 
Similarity  to  the  DEMETER  data  shown  in  Figure  3  is  readily  apparent  as  are  some  differences. 
The  most  glaring  differences  include  the  wide  empty  area  near  the  equator  (both  Figure  3  and  15 
show  a  dynamic  range  of  5  orders  of  magnitude),  the  lack  of  power  in  the  conjugate  locations, 
and  the  lack  of  any  auroral  power  spectral  density.  Part  of  the  equatorial  difference  can  be 
ascribed  to  the  failure  of  the  longitudinal  approximation  there  and  so  the  validity  of  the  Helliwell 
curves.  Another  part  is  related  to  the  second  point;  propagation  of  VLF  radiation  is  coupled  to 
magnetic  field  lines.  We  expect  to  see  power  spectral  density  near  the  conjugate  locations  of 
strikes,  but  we  also  expect  power  spectral  density  to  be  displaced  along  field  lines  (i.e. 
equatorward)  on  its  travel  from  the  ground  to  the  satellite  location.  This  effect  has  been  observed 
previously  in  DEMETER  data,  e.g.  [54],  and  will  tend  to  fill  in  the  blank  areas  near  the  equator. 
Though  trans ionospheric  absorption  was  applied,  Figure  15  still  represents  locations  within  the 
earth  ionosphere  waveguide. 
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2.7  Magnetospheric  Propagation 

As  previously  discussed  both  ducted  and  non-ducted  propagation  of  whistler  waves  in  the 
magnetosphere  have  been  observed,  or  inferred,  and  both  are  thought  to  be  important  in  different 
contexts.  Ducted  propagation  follows  field  lines  and  so,  to  first  order,  power  at  a  given  apex 
latitude  shows  up  at  an  equivalent  conjugate  latitude  with  frequency  dependent  temporal 
dispersion  and  little  damping.  Non-Ducted  propagation  is  more  complex  showing  a  strong 
coupling  between  spatial  and  spectral  dependencies.  Power  input  at  a  given  apex  latitude  can 
show  up  at  either  a  relatively  higher  or  lower  conjugate  latitude  depending  on  frequency.  In  the 
current  treatment  only  ducted  propagation  will  be  explicitly  handled,  though  it  would  certainly 
be  of  interest  to  extend  this  work  to  non-Ducted  propagation.  We  handle  it  as  follows:  field  line 
trajectories  were  calculated  at  each  comer  of  a  l°xl°  geocentric  grid  as  discussed  for  the  apex 
latitude  calculations;  pseudopower  was  calculated  within  each  grid  location;  the  pseudopower 
was  assumed  to  remain  within  the  flux  tube  defined  by  the  comers  of  each  grid  and  transported 
to  different  altitudes  within  the  apex  coordinate  system;  in  order  to  retain  per  meter  squared  units 
the  pseudopower  was  scaled  by  the  changing  cross  sectional  area  of  the  flux  tube  at  any  altitude. 
Because  we  are  interested  in  comparing  to  DEMETER  data,  an  altitude  of  660  km  will  be 
focused  upon  here  along  with  its  conjugate  location,  but  values  are  available  anywhere  along  the 
field  lines. 
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PseudoPower/km2/day,  January,  1100  hours,  5.0  -  5.1  kHz 
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PseudoPower/kryi2/dgy,  July,  2300  hours,  5.0  —  5.1  kHz 

*  660  km 


80 

60 

40 

20 

0 

-20 

-40 

-60 

-80 


d 


160  -120  -80  -40  0  40  80  120  160 

PseudoPower/km2/day,  July,  2300  hours,  5.0  —  5.1  kHz 

- -  r . —  Conjugate 


80 

60 

40 


-160  -120  -80  -40  0  40  80  120  160 

^PseudoPower/ktr^/doy,  July,  2500  hours,  5.0  -  5.1  kHz 

Sum 


B0 


20 


40 


60 


-80 


80- 
60 
40 
20 
0 

-20 
-40 
-60 
-80 


-160  -120  -80  -40  0  40  80  120  160  -160  -120  -80  -40  0  40  80  120  160 


Figure  16.  Pseudopower  Climatology  in  Low  Earth  Orbit,  Magnetospheric  Propagation 

Figure  16  shows  results  of  the  ducted  transport  calculation  applied  to  the  pseudopower 
distribution  after  Helliwell  attenuation  shown  in  Figure  15  for  January  1100  hours  local  time  and 
July  2300  hours  local  time.  Note  that  we  have  mapped  the  values  from  apex  coordinates  back  to 
a  l°xl°  geocentric  grid  at  the  altitude  in  question.  Along  with  the  model  at  660  km  (Figure  15a 
and  15b)  and  the  conjugate  location  at  660  km  (Figure  15c  and  15d)  we  also  show  an  average  of 
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the  two.  At  a  first  cut  we  think  that  the  ducted  model  is  appropriate  here  in  part  because 
DEMETER  sees  a  lot  of  power  at  5  kHz  in  the  conjugate  region  of  the  North  American  continent 
(Fig.  4d).  In  general  for  non-Ducted  propagation  we  would  expect  frequencies  below  the  LHR 
maximum  along  the  trajectory  to  undergo  magnetic  reflection  before  reaching  660  km.  As  shown 
in  Figure  1 1  and  12  there  is  significant  power  with  the  spatial  temporal  signature  of  lightning 
below  the  LHR  frequency  in  the  assumed  conjugate  region.  There  is  almost  no  lightning  below 
that  area  which  could  generate  0+  whistlers  at  the  observed  average  spectral  densities.  It  is  also 
possible  that  the  observed  radiation  is  still  propagating  in  the  longitudinal  mode  and,  although 
nonducted,  has  not  yet  become  oblique  enough  to  undergo  magnetic  reflection.  This  topic  is 
beyond  the  scope  of  the  current  investigation  and  is  subject  to  continued  study. 

Figure  15e  and  15f  show  a  marked  resemblance  to  Figure  3a  and  3d  respectively.  Outstanding 
are  the  existence  of  very  low  values  near  the  equator  in  the  night  time  summer,  which  does  not 
show  up  in  the  DEMETER  average  data,  and  that  the  data  with  intensity  below  -3 
log(pVVmTHz)  in  January  daytime  is  somewhat  broader  in  longitude  especially  near  30°  and 
135°  longitude.  While  we  have  used  -3  log(pV2/m2/Hz)  as  the  lower  limit  in  our  DEMETER  plot 
Figure  4  and  subsequent  discussion  shows  that  many  of  the  averaged  January  values  seem  to  be 
affected  by  the  instrumental  noise  floor  and/or  sensitivity.  It  is  likely  more  appropriate  to 
consider  any  value  below  -2.5  log(pV  /m  / Hz)  or  even  higher  to  be  noise  in  this  context. 

In  general  lightning  locations  are  driven  by  storm  locations,  and  while  their  occurrence  can  be 
described  climatologically,  specific  storm  and  flash  locations  are  stochastic  events.  So  even  if 
VLF  intensities  as  measured  by  DEMETER  were  entirely  due  to  lightning,  and  OTD/LIS  were 
100%  accurate  and  precise,  and  our  model  sufficiently  captured  all  features  of  the  conversion,  we 
do  not  expect  a  one  to  one  correlation  between  DEMETER  data  and  our  model.  It  is  worth  noting 
that  the  20-30  observations,  representing  800-1200  individual  spectrum,  per  grid  location  per 
month  still  only  represents  around  30  minutes  of  integrated  dwell  time.  Also  the  OTD  data  we 
are  comparing  to  have  undergone  additional  smoothing  in  order  to  show  more  continuous 
features.  For  example  the  HRMC  represents  a  2.5°  moving  average  applied  to  each  0.5°  grid, 
while  the  LRADC  involves  a  1 10  day  moving  average  centered  on  each  day  of  the  year.  We  have 
not  applied  this  level  of  smoothing  to  the  DEMETER  data.  Nevertheless,  it  seems  that  many  of 
the  gross  features  of  the  DEMETER  survey  data  can  be  explained  on  the  basis  of  this  relatively 
crude  analysis. 

2.8  Scaling  the  Model  to  DEMETER  Observations 

The  previous  analysis  has  shown  general  agreement  between  the  pseudopower  model  and  the 
averaged  DEMETER  survey  data  when  considered  in  the  appropriate  light.  We  will  now 
consider  how  best  to  scale  the  pseudopower  climatology  to  the  DEMETER  observations.  We 
know  that  not  all  of  the  spectral  density  in  the  DEMETER  data  represents  LGWs,  and  we  know 
that  many  salient  features  of  the  transfonnation  from  sferic  to  LGW  at  LEO  are  being  treated 
here  in  a  less  than  ideal  fashion.  Nevertheless,  we  will  assume  here  that  if  the  pseudopower 
climatology  is  representative  of  the  contribution  of  LGW  to  VLF  intensity  at  LEO  then  they 
should  scale  linearly.  That  is,  the  pseudopower  climatology  is  currently  on  an  arbitrary  scale,  we 
would  like  to  convert  it  to  real  units  with  a  single  scaling  factor  which  would  be  applied 
everywhere;  e.g.  pseudopower  units  *  100  =  pV  /m  / Hz.  Because  we  have  been  presenting  the 
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Scaling  Factor 


data  on  a  log  scale  this  is  equivalent  to  an  additive  scaling  in  log  units;  e.g.  log(pseudopower 
units)  +  2.  =  log(pV2/m2/Hz). 

In  determining  this  scaling  factor  we  would  like  to  only  include  observations  where  LGWs 
dominate.  Low  frequencies  where  Hiss  dominates  should  be  excluded,  as  should  frequencies 
above  the  LHR  where  quasi-electrostatic  waves  may  dominate;  from  figures  9-12  and  similar  2 
-  9  kHz  seems  reasonable  for  nighttime  and  3-7  kHz  for  daytime.  High  latitude  observations 
where  auroral  phenomena  occur  should  also  be  excluded.  As  in  Figures  10  -  12  we  choose  an 
apex  latitude  of  50°  as  a  spatial  cutoff  for  auroral  phenomena.  However  in  the  DEMETER  data 
apparently  auroral  processes  penetrate  to  lower  latitudes  at  lower  frequencies  and  during  the  day. 
These  strong  fluxes  bias  the  scaling  factor  such  that  we  have  used  a  lower  cutoff  of  4  kHz  for 
nighttime  and  5  kHz  for  daytime.  This  choice  represents  a  compromise  between  excluding  all 
auroral  phenomena  and  including  regions  where  we  expect  lightning  to  be  important.  In  addition 
we  would  like  to  be  cognizant  of  known  weaknesses  in  our  modeling  approach.  The  longitudinal 
approximation  invoked  in  the  Helliwell  calculations  breaks  down  for  small  dip  angles  and  so  we 
exclude  apex  latitudes  below  20°  from  our  determination.  Observational  limitations  such  as  the 
sampling  density  relative  to  the  signal  variability  and  the  effects  of  an  instrumental  cutoff  should 
also  be  kept  in  mind.  Finally,  the  correlation  will  be  computed  with  DEMETER  data  averaged 
over  a  2.5°  spatial  resolution,  in  order  to  smooth  over  some  of  the  stochastic  spatial  variation  in 
strikes  and  to  promote  overall  similarity  to  the  smoothing  applied  to  the  optical  flash  data  set. 

Within  the  set  of  data  limited  as  described  above  we  calculate  a  scaling  factor  from  the  fraction 
of  data  representing  relative  rank  0.6  -  0.9  by  spectral  density.  This  last  will  help  to  exclude 
outliers  and  minimize  the  contribution  of  the  instrumental  cutoff,  while  emphasizing  locations 
where  there  are  many  0+  whistlers  (as  seen  in  Figure  4).  Because  we  are  considering  the 
correlation  between  log  scaled  data  the  importance  of  the  higher  ranks  are  magnified.  Figure  17a 
shows  the  scaling  factors  for  each  month  of  the  year  at  selected  frequencies,  day  and  night. 
Figure  17b  shows  the  associated  correlation  coefficients  between  the  data  and  model. 


Figure  17.  Climatology  Scaling  Factors  and  Correlation  Coefficients 

Figure  17a  reveals  a  systematic  temporal  variation  in  the  scaling  factor,  with  a  large  peak  in 
Northern  summer  and  a  smaller  peak  in  Southern  summer.  This  is  due  at  least  in  part  to  the 
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overly  simplified  transionospheric  attenuation  treatment  where  only  the  solar  zenith  angle  as  a 
function  of  time  defines  a  season.  It  is  well  established  that  the  ionosphere  has  a  distinct  seasonal 
cycle  which  we  are  ignoring  here.  The  same  monthly  trend  is  evident  day  or  night.  It  is  also  clear 
that  the  daytime  and  nighttime  values  fall  within  different  distributions,  and  the  scaling  factor  is 
roughly  an  order  of  magnitude  higher  during  the  day.  That  is,  given  only  a  single  scaling  factor 
based  on  nighttime  data,  the  model  would  underestimate  the  spectral  densities  by  a  systematic 
factor  during  the  day.  This  most  likely  relates  to  our  simple  transionospheric  transfer  function, 
but  could  also  have  a  contribution  from  the  OTD/LIS  data  that  our  model  originated  from,  which 
might  show  different  sensitivities  between  day  and  night.  The  correlation  coefficients  shown  in 
Figure  17b  reveal  a  very  similar  temporal  variation.  This  suggests  another  possible  contributing 
factor  to  the  monthly  trend:  there  is  increased  lightning  activity  at  high  latitudes  during  summer, 
while  more  lightning  is  equatorial  near  the  equinox  (driven  by  the  solar  zenith  angle  differences 
between  the  seasons).  Lightning  at  high  latitudes  is  converted  to  LGW  more  efficiently  (due  to 
higher  magnetic  dip  angle,  smaller  electron  column  densities,  and  possibly  additional  ducts), 
which  provides  a  larger  signal  above  background  and  mitigates  some  of  the  stochastic 
component  in  the  signal.  Overall,  the  correlation  coefficients  shown  in  Figure  17b  reveal  values 
consistent  with  a  significant  relationship  between  the  model  and  the  data,  while  still  not 
accounting  for  all  of  the  variance. 

It  would  have  been  preferable  to  have  a  single  scaling  factor  between  the  pseudopower 
representing  our  model  and  the  DEMETER  spectral  density,  but  Figure  17a  shows  that  it  is 
inadvisable.  This  is  not  unexpected  given  the  first  cut  treatment  given  here,  but  may  still  be 
achievable  in  the  future.  At  minimum  different  factors  should  be  used  for  daytime  and  nighttime, 
but  given  the  apparent  systematic  nature  of  the  monthly  trends  in  both  the  scaling  factor  and 
correlation  coefficient  across  frequencies  and  local  time  it  seems  reasonable  to  use  a  different 
scaling  factor  for  each  month,  day  and  night.  The  final  values  are  averaged  over  relevant 
frequencies  and  shown  in  Table  1. 


Table  1.  Frequency  Averaged  Scaling  Factors 


Jan. 

Feb. 

Mar. 

Apr. 

May 

June 

July 

Aug. 

Sep. 

Oct. 

Nov. 

Dec. 

Mean 

Day 

4.30 

4.44 

4.32 

4.35 

4.52 

4.65 

4.83 

4.81 

4.70 

4.40 

4.35 

4.28 

4.50 

Night 

3.75 

3.76 

3.67 

3.61 

3.76 

3.86 

3.97 

3.94 

3.84 

3.69 

3.58 

3.71 

3.76 

4.  DISCUSSION 

Figure  18  gives  an  overview  of  the  model  values  compared  to  the  DEMETER  averages  for  5 
kHz,  day  and  night,  January,  April,  July,  and  October.  The  equatorial  region  is  clearly  in  need  of 
a  better  treatment  of  transionospheric  attenuation  and/or  sferic  to  LGW  transfonnation, 
especially  at  night.  This  fact  is  consistent  with  the  breakdown  of  the  longitudinal  approximation 
there,  but  also  likely  represents  the  need  for  a  seasonally  varying  ionosphere  and  perhaps 
adjustment  for  the  known  heterogeneity  of  the  equatorial  ionosphere  (e.g.  scintillation  and 
plasma  depletions).  In  the  recent  reanalysis  of  the  Helliwell  curves  [60]  the  only  glaring 
difference  between  the  original  curves  and  the  recalculated  values  based  on  the  IRI  and  a  full 
wave  model  is  a  flattening  out  of  the  absorption  for  nighttime  values  below  20°  latitude  (Figure  5 
in  that  work),  and  such  an  adjustment  would  certainly  improve  the  agreement  but  still  not  close 
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the  gap  entirely.  The  need  for  an  improved  earth-ionosphere  waveguide  can  perhaps  be  seen  in 
the  regions  away  from  active  lightning.  As  one  moves  longitudinally  away  from  an  active 
lightning  region  along  a  specific  apex  latitude  the  DEMETER  data  seems  to  drop  off  more 
quickly  than  the  model.  This  could  be  driven  by  absorption  in  the  waveguide,  or  alternatively  by 
the  increasing  obliquity  of  the  wave  with  respect  to  the  ionospheric  plasma  gradient  with 
increasing  distance  from  the  stroke.  In  addition  our  simplified  treatment  implicitly  assumes  that 
the  spectrum  of  the  sferic  does  not  change  with  distance,  which  is  not  true  at  least  within  the 
waveguide  as  discussed  previously.  An  initial  look  for  regions  where  this  effect  might  be 
important  (e.g.  over  the  Pacific  Ocean)  surprisingly  revealed  no  large  differences  in  spectrum. 
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Figure  18.  Spatial  Plots  of  the  Scaled  Climatology  and  DEMETER  Averages 


The  stochastic  nature  of  lightning  and  storms  relative  to  satellite  observations  is  also  evident  in 
Figure  18.  Even  within  regions  of  high  lightning  activity  the  model  is  smooth.  The  values 
originate  in  a  highly  smoothed  OTD  data  set  and  become  even  smoother  due  to  our  inverse 
distance  treatment.  While  a  more  detailed  waveguide  treatment  may  improve  the  model  it  will 
still  introduce  a  similar  smoothing  effect.  In  order  to  capture  the  heterogeneity  of  the  DEMETER 
data  it  would  be  necessary  to  treat  individual  storms  as  events  lasting  a  finite  time  and  lightning 
as  a  stochastic  process  where  an  individual  stroke  was  drawn  from  a  representative  power 
distribution.  In  addition  the  conversion  from  waveguide  power  to  LGW  power  may  have  a 
stochastic,  or  at  least  a  spatial  and  temporal,  component  as  it  may  depend  on  local  variations  in 
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the  ionosphere  (e.g.  horizontal  gradients  in  electron  density),  geomagnetic  dip  angle,  or  distance 
traveled  in  the  waveguide.  Even  if  these  approaches  were  implemented  it  would  still  be 
necessary  to  model  the  satellite  as  sampling  this  distribution  as  it  orbits.  That  is,  even  with  a 
mission  like  DEMETER  where  survey  mode  data  was  captured  almost  continuously  for  five 
years,  each  l°xl°  location  was  only  sampled  for  roughly  30  minutes  in  any  given  month.  There 
was  a  large  range  in  the  values  over  a  specific  location  for  a  given  overpass,  and  this  can  be  seen 
in  the  almost  vertical  striations  observed  in  the  DEMETER  plots.  These  variations  must  depend 
in  part  on  whether  there  was  active  stonn  below  the  satellite  or  at  some  distance,  and  when 
individual  strikes  were  integrated  into  the  survey  mode  average  spectra.  We  are  implicitly 
assuming  that  the  true  distribution  was  sampled  sufficiently  but  it  would  be  of  interest  to  test  that 
assumption. 

Our  climatology  also  extends  to  local  times  not  sampled  by  the  DEMETER  mission.  Figure  19 
shows  the  spectral  density  at  5  kHz  plotted  vs.  apex  radius  for  some  of  our  local  time  bins  with 
color  representing  different  months  of  the  year  where  violet  -t>  red  corresponds  to  January  -» 
December.  It  is  interesting  that  Figure  19f,  which  overlaps  with  the  nighttime  DEMETER  orbit, 
also  has  the  most  spectral  density.  This  reflects  a  combination  of  flash  density  and  ionospheric 
attenuation.  Typically  the  peak  spectral  density  was  in  the  2000  -  2200  hours  local  time  bin  for  a 
given  apex  radius,  frequency,  and  month.  This  is  expected  as  it  contains  the  peak  flash  rate 
(Figure  7)  under  a  nighttime  ionosphere.  Figure  19  makes  it  evident  that  spectral  density  from 
lightning  is  a  complicated  function  of  the  apex  radius,  month,  and  local  time  as  different  land 
masses  experience  different  meteorological  and  ionospheric  conditions. 
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Figure  19.  Spectral  Density  of  the  Scaled  Climatology  vs.  Apex  Radius 


5.  CONCLUSIONS 

Starting  from  OTD/LIS  lightning  data  representing  1995  -  2005,  we  have  constructed  a 
climatology  of  VLF  spectral  density  from  lightning  at  LEO.  Key  assumptions  include:  1)  a  linear 
relationship  between  optical  flash  rate  and  VLF  power  flux;  2)  VLF  attenuation  goes  as  the 
inverse  distance  in  the  earth-ionosphere  waveguide  (includes  waveguide  attenuation  and  sferic- 
whistler  coupling);  3)  transionospheric  attenuation  can  be  reasonably  represented  by  the 
Helliwell  curves;  4)  0+  whistler  energy  travels  along  field  lines  and  shows  up  in  the  conjugate 
hemisphere  without  attenuation;  and  5)  the  5  years  of  DEMETER  survey  data  reasonably 
represents  the  spectral  density  in  LEO.  Values  are  available  for  each  month  of  the  year  and  every 
two  hours  of  local  time  on  a  l°xl°  spatial  grid.  The  so  constructed  climatology  shows 
complicated  variation  as  a  function  of  time  and  space. 

This  climatology  represents  a  baseline  treatment  and  there  are  some  obvious  improvements  that 
could  be  tested  and  implemented.  In  addition  it  should  be  kept  in  mind  that  the  observations  to 
which  this  model  is  tied  were  mostly  taken  during  the  occurrence  of  an  extreme  solar  minimum. 
There  are  parameters  which  could  affect  observed  spectral  power  densities  and  which  would  also 
be  affected  by  solar  minimum  conditions;  such  as  the  phase  space  density  of  energetic  electrons 
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throughout  the  inner  magnetosphere,  and  the  location  of  the  plasmapause.  We  now  provide  a 
brief  summary  of  the  models  used  for  each  organizational  area  of  the  framework  including 
weaknesses  and  potential  for  improvement. 

5.1  Terrestrial  VLF  Source  Model 

The  Terrestrial  VLF  source  model  sprang  from  a  combination  of  widely  available  LIS/OTD  flash 
rate  data  sets  provided  by  the  Global  Hydrology  and  Climate  Center's  Lightning  Team.  This 
same  archive  also  contains  time  series  data  which  could  be  used  to  infer  the  year  to  year  spatial 
and  temporal  variations  associated  with  the  averages.  Conversion  to  power  flux  was  done  via  the 
assumption  that  flash  rate  is  directly  proportional  to  VLF  power  emitted.  Spectral  density  was 
then  derived  by  assuming  an  average  spectral  content  per  flash.  The  primary  drawback  of  this 
procedure  is  that  each  flash  in  every  stonn  is  considered  equivalent.  It  would  be  better  to  use  a 
data  set  that  included  infonnation  about  the  distribution  of  power  among  strokes  spatially, 
temporally,  and  even  spectrally.  Ground  based  lightning  networks  are  currently  operating  which 
give  global  information  on  the  spatial  and  temporal  variation  of  VLF  power,  but  they  are  still 
mostly  in  the  research  and  development  phase.  These  same  networks  can  give  spectral 
information  as  a  function  of  distance  from  the  strike.  It  would  be  of  interest  to  use  data  from  such 
networks  to  at  least  test  the  assumptions  listed  above. 

5.2  Earth-Ionosphere  Waveguide  Propagation  Model 

The  Earth-Ionosphere  waveguide  propagation  model  employed  here  is  the  simplest  form,  an 
inverse  square  law.  Within  line  of  sight  viewing  of  the  lightning  strike  from  the  bottom  of  the 
ionosphere  this  is  likely  a  good  approximation.  Simple  geometry  tells  us  that  from  an  altitude  of 
100  km  the  line  of  sight  viewing  is  approximately  cos'^OSSO/hdSO)  =  10°  or  about  1000  km. 

This  is  very  similar  to  the  maximum  distance  for  sferic  to  whistler  conversion  estimated  from  in 
situ  rocket  and  satellite  observations  [32,  33,  51].  The  main  drawback  to  this  approach  is  that 
there  is  no  frequency  dependent  attenuation.  It  is  likely  that  our  approach  underestimates 
attenuation  at  intennediate  distances  for  frequencies  below  about  5  kHz  and  overestimates  it 
above  (see  the  mode  theory  discussion  in  Section  6  and  [56])  but  some  of  this  effect  is  mitigated 
by  the  much  larger  ionospheric  absorption  at  higher  frequencies  and  the  increasing  reflection 
coefficient  with  increasing  wave  obliquity  with  respect  to  the  lower  ionospheric  boundary.  We 
most  likely  underestimate  attenuation  at  long  distances  for  all  frequencies. 

Short  of  a  full  waveguide  mode  treatment  of  the  propagation  it  should  be  possible  to  transition 
from  inverse  square  attenuation  to  a  frequency  dependent  dB/Mm  model  at  a  distance  near  1000 
km.  Different  values  should  be  used  for  day  and  night  due  the  large  difference  in  the  lower 
ionospheric  boundary  altitude  seen  by  VLF  radiation.  It  is  also  probably  true  that  the  sferic  to 
whistler  conversion  becomes  less  efficient  at  large  distances  as  the  wave  vector  becomes  mostly 
horizontal  to  the  earth's  surface.  Finally,  there  should  be  a  geomagnetic  intensity  and  dip  angle 
dependence.  These  effects  could  be  incorporated  with  a  spatial  dependence  to  the  dB/Mm 
treatment.  Spectrally  dependent  waveguide  transport  may  be  responsible  for  the  fact  that  the 
overall  spectral  peak  appears  near  4  kHz.  The  peak  might  be  expected  at  a  lower  frequency  based 
on  the  incipient  sferic  spectrum  and  the  steep  rise  in  absorption  as  frequency  decreases.  An 
ameliorating  factor  might  be  the  smaller  attenuation  experienced  by  higher  frequencies  in  the 
waveguide.  In  addition  there  is  very  strong  absorption  below  about  2  kHz  due  to  the  so  called 
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cutoff  frequency  of  the  earth-ionosphere  waveguide,  e.g.  [62].  This  cutoff  effect  could  drive  the 
minimum  seen  at  1.6  kHz  in  Figure  9a  which  shows  the  seasonality  of  lightning.  The  other 
minimum  below  1  kHz  is  most  likely  related  to  the  proton  cyclotron  frequency  in  the  ionosphere 
below  which  whistlers  are  said  not  to  propagate.  The  proton  cyclotron  frequency  is  proportional 
to  the  magnitude  of  the  magnetic  field  and  so  varies  spatially,  but  is  typically  around  600  Hz  in 
regions  of  interest.  Neither  of  these  effects  are  currently  being  treated  and  so  power  at 
frequencies  below  2  kHz  is  probably  overestimated  in  our  climatology. 

5.3  Trans-Ionospheric  Attenuation  Model 

The  Trans-Ionospheric  attenuation  model  employed  here  originates  in  work  by  Helliwell  [40].  In 
that  work  two  simplified  electron  profiles  in  the  ionosphere  were  assumed,  one  for  daytime  and 
another  for  nighttime  conditions.  The  calculation  proceeded  in  two  parts  one  for  the  lower 
ionosphere  (60  -  200  km)  where  electron  neutral  collisions  dominate  and  the  second  for  the  F 
region  (200  -  1500  km)  where  Coulomb  interactions  between  electrons  and  ions  dominate.  A 
single  vertical  profile  for  each  of  total  collisions  and  Coulomb  collisions  was  employed. 
Absorption  was  calculated  for  vertically  incident  waves  along  the  wave  normal  direction  using  a 
refractive  index  representing  the  Quasi-Longitudinal  (QL)  approximation.  As  previously 
discussed  the  major  deficiency  is  the  lack  of  known  systematic  variation  in  ionospheric  profiles. 
The  calculation  is  represented  in  some  detail  in  [40]  and  was  recently  reproduced  omitting  the 
QL  assumption,  and  then  again  with  a  full  wave  code  [60],  which  showed  the  overall  validity  of 
the  Helliwell  approach.  When  compared  to  average  ionospheric  conditions  calculated  from  the 
IRI  the  major  difference  was  due  to  a  sharper  cutoff  in  electron  densities  below  80  km  seen  in  the 
IRI  nighttime  profile,  this  difference  lead  to  somewhat  lower  absorption.  It  would  certainly  be  of 
interest  to  use  IRI  based  ionospheric  profiles,  of  whatever  fidelity  is  computationally  feasible.  In 
leu  of  perfonning  the  large  number  of  full  wave  calculations  necessary  to  cover  the  range  of 
dimensions  addressed  here  (longitude,  latitude,  time  of  day,  frequency  =  360x120x12x205)  it 
might  be  expeditious  to  use  the  calculation  technique  exactly  as  outlined  by  Helliwell,  but  with 
time  and  space  electron  profiles  from  the  IRI;  especially  in  light  of  the  close  agreement  between 
it  and  the  full  wave  calculation  shown  in  [60]  (Figure  2  of  that  work)  when  identical  input 
parameters  are  utilized. 

5.4  Magnetospheric  Propagation  Model 

The  Magnetospheric  propagation  model  employed  here  is  pure  ducted.  That  is,  we  assume  field 
aligned  propagation  only,  while  in-situ  magnetospheric  observations  indicate  that  most  whistlers 
are  likely  nonducted,  e.g.  [63].  It  is  known  that  whistler  waves  propagating  in  non-ducted  mode 
through  the  ionosphere  maintain  vertical  wave  nonnal  directions  due  to  the  strong  gradient  in 
electron  density;  but  their  group  velocity  direction  is  not  vertical  in  this  region  due  to  the 
influence  of  the  geomagnetic  field.  We  have  modeled  transport  through  the  ionosphere  up  to 
LEO  as  ducted  (field  aligned)  propagation  while  ducts  are  thought  to  remain  coherent  only  down 
to  1000  km  or  so.  Thus  transport  of  0+  whistler  up  to  LEO  altitudes  should  likely  be  non-ducted 
and  arrive  somewhere  between  radial  and  field  aligned  with  a  frequency  dependent  component. 
This  picture  is  borne  out  by  observational  studies  as  previously  described  [54].  Given  the  fidelity 
of  this  initial  study  such  a  correction  is  not  expected  to  be  important,  but  could  be  simply 
implemented  in  the  future. 
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DEMETER  observations  show  large  amounts  of  power  with  the  spatial/temporal  character  of 
lightning  at  and  below  5  kHz  in  regions  where  there  is  very  little  lightning  activity,  but  where 
there  is  activity  in  the  conjugate  region.  Non-ducted  transport  calculations  at  these  frequencies 
indicate  that  such  waves  should  undergo  MR  above  DEMETER  orbit  altitudes.  So  it  seems  that 
ducted  transport  is  important.  This  may  include  one  sided  ducting  associated  with  one  way 
plasma  gradients  such  as  the  plasmapause  [64]  as  well  as  the  traditional  field  aligned 
perturbations.  DEMETER  observations  also  show  LHR  hiss  in  these  same  conjugate  regions, 
which  is  undoubtedly  associated  with  non-ducted  waves.  So,  non-ducted  transport  is  also 
important.  The  LHR  associated  waves  could  be  high  enough  frequency  to  not  have  undergone 
MR  before  reaching  LEO.  The  one  axis  electric  field  measurements  considered  here  cannot 
distinguish  between  longitudinal  and  oblique  propagation.  There  is  some  suggestion  in  the 
literature  that  the  chance  of  ducted  propagation  has  a  latitudinal  dependence,  with  predominantly 
non-ducted  transport  being  available  below  L=1.6  or  so  [65,  66].  Comparing  the  model  and 
observations  for  July  nighttime  at  5  kHz  shown  in  Figure  18  we  see  that  the  low  latitude  leg  of 
the  North  American  continent  (Mexico  -  Panama)  shows  up  in  the  conjugate  model  intensities 
but  not  in  the  observations;  consistent  with  this  hypothesis.  It  seems  likely  that  the  LHR  hiss 
observed  in  conjugate  regions  has  its  origins  in  non-ducted  waves  originating  from  lower 
latitudes,  while  the  spectral  density  around  5  kHz  has  its  origins  in  ducted  waves  from  the 
conjugate  field  line  location. 

In  practice  the  assumption  of  ducted  transport  simplifies  the  calculations  greatly.  This  is  because 
there  is  no  spatially  dependent  frequency  dispersion.  While  in  non-ducted  transport  the  wave 
embarks  on  a  whole  new  spatial  trajectory  at  each  MR  until  it  settles  at  an  apex  radius 
characteristic  of  its  frequency  (which  has  an  equatorial  LHR  approximately  equal  to  the  wave 
frequency  according  to  [67]).  Thus  for  ducted  transport  a  given  location  in  the  apex  coordinate 
system  only  has  contributions  from  other  locations  along  that  field  line,  but  for  nonducted 
transport  a  given  location  can  have  contributions  from  a  large  number  of  other  locations.  In  this 
initial  treatment,  and  where  we  are  only  trying  to  represent  spectral  density  at  LEO,  we  have 
found  a  pure  ducted  treatment  expedient  and  sufficient.  However,  in  order  to  fill  the  inner 
magnetosphere  with  waves  from  lightning  it  would  be  necessary  to  make  some  assumption  about 
at  least  the  relative  occurrence  of  ducted  and  non-ducted  waves.  As  in  Able  and  Thorne  [1998], 
where  it  was  assumed  that  each  accounted  for  roughly  the  same  amount  of  spectral  density  and 
the  same  relative  occurrence.  While  such  an  estimate  is  beyond  the  scope  of  this  initial 
investigation  we  will  note  that  the  overall  spatial  temporal  spectral  infonnation  within  the 
DEMETER  data  set  could  give  some  insight  in  to  the  relative  proportions  for  reasons  already 
outlined  above. 

The  ultimate  purpose  of  this  investigation  was  to  give  estimates  of  the  VLF  power  in  LEO  at  a 
two  hour  cadence  for  each  month  of  the  year.  However,  in  order  to  characterize  the  sources  and 
sinks  of  energetic  electrons  in  the  radiation  belts  the  spectral  power  densities  derived  here  must 
be  used  as  input  to  a  calculation  of  pitch  angle  diffusion  coefficients  which  could  then  be  used  to 
estimate  particle  lifetimes.  Although  we  have  used  only  ducted  transport  calculations  here,  in 
practice  these  spectral  power  densities  could  be  used  as  input  to  either  ducted  or  nonducted 
transport  calculations  since  the  scalings  and  correlations  are  based  on  total  power.  The 
climatologies  reported  here  represent  a  step  forward  from  the  single  numbers  used  in  Able  and 
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Thorne  [10],  or  even  the  L  and  AE  dependent  values  derived  by  Meredith  et  al.  [23].  Work  on 
converting  these  spectral  power  density  estimates  into  radiation  belt  particle  lifetimes  is  ongoing. 
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